package analysis;


/**
 * @author Mathieu Rousseau
 * @email mathieu.rousseau3@mail.mcgill.ca
 */


import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.text.DateFormat;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Date;

import utilities.PDButility;


/** ------------------------------------------------------------------------------------------------------------------
 * start of Microcosm class
 * ------------------------------------------------------------------------------------------------------------------  **/
public class Microcosm {


	/*  global variables  */
	private final String inputPath = "Microcosm/Input/";
	private final DateFormat dateFormatOutputFolder = new SimpleDateFormat("yyyy-MM-dd-HH-mm-ss");
	private final String outputFolderDate = dateFormatOutputFolder.format(new Date());
	private String outputFolder = "Microcosm/Output/" + outputFolderDate;
	private String outputPath = "Microcosm/Output/";
	private double radius = 0.1;
	private double radiusSquared = radius*radius;
	private int resolution = 1;
	private PDButility pdb = new PDButility();
	private int numberOfStructures = 0;
	private int numberOfPoints = 0;
	private String experimentName = "none";
	private String fragmentsName = "default";
	private ArrayList<double[][]> structures;
	private int[] fragments;
	private int minBase = -1;
	private int maxBase = -1;
	private int numberOfBases = -1;
	private ArrayList<double[]> baseDensities;
	double averageBaseDensity[];
	double stdevBaseDensity[];
	private ArrayList<double[]> baseCondensations;
	double averageBaseCondensation[];
	double stdevBaseCondensation[];
	private ArrayList<double[]> baseLoopings;
	double averageBaseLooping[];
	double stdevBaseLooping[];
	/*  end of global variables  */
	
	
	/*  start of getStructuresFromFile method  */
	private void getStructuresFromFile(){
		//  declare and initialize the list that will hold all of the structures
		structures = new ArrayList<double[][]>(numberOfStructures);
		for(int i=0; i<numberOfStructures; i++){
			double structureInitializer[][] = new double[numberOfPoints][numberOfPoints];
			structures.add(structureInitializer);
		}

		
		//  define the filename of each of the structures to be read from file
		String structureFilename;
		//  use the convention that the structure file names start at index 1
		String suffix = "";
		if(!experimentName.equals("none")){
				suffix += "_" + experimentName;
		}
		for(int i=1; i<=numberOfStructures; i++){
			if(i <= 9){
				structureFilename = "structure_000" + i + suffix + ".pdb";
			}
			else if(i <= 99){
				structureFilename = "structure_00" + i + suffix + ".pdb";
			}
			else if(i <= 999){
				structureFilename = "structure_0" + i + suffix + ".pdb";
			}
			else{
				structureFilename = "structure_" + i + suffix + ".pdb";
			}
//			System.out.println("Tyring to read structure file: " + inputPath + structureFilename);
			
			//  get the structure from file and store it in the list of structures
			//  remember that the indices in the array list start at zero (0) and not one (1) like the filenames
			structures.set(i-1, pdb.getStructure(inputPath + structureFilename));
		}
	}/*  end of getStructuresFromFile method  */
	
	
	/*  start of getFragmentsFromFile method  */
	private void getFragmentsFromFile(){
		fragments = new int[numberOfPoints];
		
		//	open and read the fragments data file
		String line = null;
		try{
			BufferedReader br = new BufferedReader(new FileReader(inputPath + fragmentsName));			
			int lineCounter = 0;
			while((line = br.readLine()) != null){
				fragments[lineCounter] = Integer.parseInt(line);
				lineCounter++;
			}
			if(lineCounter != numberOfPoints){
				System.out.println("Error: the number of fragments (" + lineCounter + ") does NOT equal the number of points (" + numberOfPoints + ").");
				System.out.println("System.exit(0)");
				System.exit(0);
			}
//			else{
//				System.out.println("The number of fragments (" + lineCounter + ") does equal the number of points (" + numberOfPoints + ").");
//			}
			br.close();
		}catch (IOException e) {
			e.printStackTrace();
		}
		
		
		//  make sure that the base indices are monotonically increasing
		for(int i=0; i<numberOfPoints-1; i++){
			if(!(fragments[i] < fragments[i+1])){
				System.out.println("Error: the base indices are NOT monotonically increasing (" + fragments[i] + " is NOT less than " + fragments[i+1] + ").");
				System.out.println("System.exit(0)");
				System.exit(0);
			}
		}
		
		
		//  determine the number of bases in one structure
		minBase = Integer.MAX_VALUE;
		maxBase = Integer.MIN_VALUE;
		for(int i=0; i<fragments.length; i++){
			int currentValue = fragments[i];
			if(currentValue < minBase){
				minBase = currentValue;
			}
			if(currentValue > maxBase){
				maxBase = currentValue;
			}
		}
		numberOfBases = maxBase - minBase + 1;
//		System.out.println("\tminBase = " + minBase + "\n\tmaxBase = " + maxBase + "\n\tNumber of Bases = " + numberOfBases);
		
		
		//  Display the fragments read from the file
//		for(int i=0; i<fragments.length; i++){
//			System.out.println("[" + (i+1) + "]\t" + fragments[i]);
//		}
	}/*  end of getFragmentsFromFile method  */
	
	
	/*  start of getBaseCoordinates method  */
	private double[] getBaseCoordinates(int structureIndex, int baseIndex){
		double[] baseCoordinates = new double[3];
		boolean alreadySetBaseCoordinates = false;
		
		
		//  find the endpoints of the fragment containing the base
		boolean foundFragmentContainingBase = false;
		int startPointBaseIndex = -1;
		int endPointBaseIndex = -1;
		int fragmentIndex = -1;
		for(int i=0; i<numberOfPoints-1; i++){
			if((fragments[i] <= baseIndex) && (baseIndex < fragments[i+1])){
				fragmentIndex = i;
				startPointBaseIndex = fragments[fragmentIndex];
				endPointBaseIndex = fragments[fragmentIndex+1];
				foundFragmentContainingBase = true;
				break;
			}
		}
		//  if cannot find the base then investigate
		if(!foundFragmentContainingBase){
			//  in the case that the base corresponds exactly to the endpoint of the last fragment, set the coordinates automatically
			if(baseIndex == maxBase){
//				System.out.println("\tLooking for the coordinates of the base (" + baseIndex + ") corresponding to the endpoint of the last fragment (a.k.a. maxBase (" + maxBase + ")).");
				fragmentIndex = numberOfPoints-1;
				startPointBaseIndex = fragments[fragmentIndex];
				endPointBaseIndex = fragments[fragmentIndex];
				foundFragmentContainingBase = true;
				baseCoordinates[0] = structures.get(structureIndex)[fragmentIndex][0];
				baseCoordinates[1] = structures.get(structureIndex)[fragmentIndex][1];
				baseCoordinates[2] = structures.get(structureIndex)[fragmentIndex][2];
				alreadySetBaseCoordinates = true;
			}
			else{
				System.out.println("Error: could not find the fragment containing base " + baseIndex + "\nSystem.exit(0)");
				System.exit(0);
			}
		}
		else{
//			System.out.println("\tbase " + baseIndex + " is located in fragment " + (fragmentIndex+1));
//			System.out.println("\t\tstart point " + startPointBaseIndex + " has coordinates\tx:" + structures.get(structureIndex)[fragmentIndex][0] + 
//					"\ty:" + structures.get(structureIndex)[fragmentIndex][1] + "\tz:" + structures.get(structureIndex)[fragmentIndex][2]);
//			System.out.println("\t\tend point " + endPointBaseIndex + " has coordinates\tx:" + structures.get(structureIndex)[fragmentIndex+1][0] + 
//					"\ty:" + structures.get(structureIndex)[fragmentIndex+1][1] + "\tz:" + structures.get(structureIndex)[fragmentIndex+1][2]);
		}
		
		
		//  calculate the coordinates of the base
		if(!alreadySetBaseCoordinates){
			double ratio = (double)(baseIndex - startPointBaseIndex) / (double)(endPointBaseIndex - startPointBaseIndex);
//			System.out.println("\tThe ratio = " + ratio);
			//  if the ratio is equal to zero then the base coordinates correspond exactly to the coordinates of the start point of the fragment
			if(ratio == 0){
				baseCoordinates[0] = structures.get(structureIndex)[fragmentIndex][0];
				baseCoordinates[1] = structures.get(structureIndex)[fragmentIndex][1];
				baseCoordinates[2] = structures.get(structureIndex)[fragmentIndex][2];
			}
			//  if the ratio is equal to one then the base coordinates correspond exactly to the coordinates of the end point of the fragment (which is exaclty the start point of the following fragment)
			else if(ratio == 1){
				baseCoordinates[0] = structures.get(structureIndex)[fragmentIndex+1][0];
				baseCoordinates[1] = structures.get(structureIndex)[fragmentIndex+1][1];
				baseCoordinates[2] = structures.get(structureIndex)[fragmentIndex+1][2];
			}
			//  if the coordinates have not already been set, then need to set the coordinates of the base using the value of the ratio
			else{
				baseCoordinates[0] = structures.get(structureIndex)[fragmentIndex][0] + (ratio * (structures.get(structureIndex)[fragmentIndex+1][0] - structures.get(structureIndex)[fragmentIndex][0]));
				baseCoordinates[1] = structures.get(structureIndex)[fragmentIndex][1] + (ratio * (structures.get(structureIndex)[fragmentIndex+1][1] - structures.get(structureIndex)[fragmentIndex][1]));
				baseCoordinates[2] = structures.get(structureIndex)[fragmentIndex][2] + (ratio * (structures.get(structureIndex)[fragmentIndex+1][2] - structures.get(structureIndex)[fragmentIndex][2]));
			}
		}
		
		
		//  display the coordinates for the base
//		System.out.println("\tThe coordinates for base " + baseIndex + " are:");
//		System.out.println("\t\tx:\t" + baseCoordinates[0]);
//		System.out.println("\t\ty:\t" + baseCoordinates[1]);
//		System.out.println("\t\tz:\t" + baseCoordinates[2]);
		
		
		return baseCoordinates;
	}/*  end of getBaseCoordinates method  */
	
	
	/*  start of getDistanceSquaredBetweenBases method  */
	private double getDistanceSquaredBetweenBases(int structureIndex, int firstBaseIndex, int secondBaseIndex){
		double[] firstBaseCoordinates = getBaseCoordinates(structureIndex, firstBaseIndex);
		double[] secondBaseCoordinates = getBaseCoordinates(structureIndex, secondBaseIndex);
		
		
		double distanceSquared = 0;
		//  first for the x-coordinates
		double term = secondBaseCoordinates[0] - firstBaseCoordinates[0];
		term *= term;
		distanceSquared += term;
		//  next for the y-coordinates
		term = secondBaseCoordinates[1] - firstBaseCoordinates[1];
		term *= term;
		distanceSquared += term;
		//  last for the z-coordinates
		term = secondBaseCoordinates[2] - firstBaseCoordinates[2];
		term *= term;
		distanceSquared += term;
//		System.out.println("The distance between the two bases (" + (firstBaseIndex+minBase) + ", " + (secondBaseIndex+minBase) +
//				" is " + Math.sqrt(distanceSquared));
		
		//  make sure that the distance is strictly greater than zero
		if(distanceSquared <= 0){
			System.out.println("Error: the distance between the two bases (" + firstBaseIndex + ", " + secondBaseIndex +
					") is equal to zero (0).\nSystem.exit(0)");
			System.exit(0);
		}
		
		return distanceSquared;
	}/*  end of getDistanceSquaredBetweenBases method  */

	
	/*  start of getDistanceBetweenBases method  */
	private double getDistanceBetweenBases(int structureIndex, int firstBaseIndex, int secondBaseIndex){
		return Math.sqrt(getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex));
	}/*  end of getDistanceBetweenBases method  */
	
	
	/*  start of calculateAverages method  */
	private double[] calculateAverages(ArrayList<double[]> arrays){
		//  calculate the average
		double[] averages = new double[0];
		if(arrays.size() == 0){
			System.out.println("Warning: trying to calculate the average from an empty set of arrays.");
		}
		else{
			int sizeOfOneArray = arrays.get(0).length;		//  know that there exists at least one element in 'arrays' object
			averages = new double[sizeOfOneArray];
			for(double[] currentArray: arrays){
				for(int i=0; i<sizeOfOneArray; i++){
					averages[i] += currentArray[i];
				}
			}
			for(int i=0; i<sizeOfOneArray; i=i+resolution){
				averages[i] /= (double)arrays.size();		//  know that there exists at least one element in 'arrays' object
			}			
		}
		return averages;
	}/*  end of calculateAverages method  */
	
	
	/*  start of calculateStandardDeviations method  */
	private double[] calculateStandardDeviations(ArrayList<double[]> arrays, double[] averages){
		//  calculate the standard deviation
		double[] stdevs = new double[0];
		if(arrays.size() < 2){
			System.out.println("Warning: trying to calculate the standard deviation from a set of arrays with less than two elements.");
		}
		else{
			int sizeOfOneArray = arrays.get(0).length;		//  know that there exists at least one element in 'arrays' object
			stdevs = new double[sizeOfOneArray];
			for(double[] currentArray: arrays){
				for(int i=0; i<sizeOfOneArray; i++){
					double term = currentArray[i] - averages[i];
					term *= term;
					stdevs[i] += term;
				}
			}
			for(int i=0; i<sizeOfOneArray; i=i+resolution){
				stdevs[i] /= (double)(arrays.size()-1);		//  know that there exists at least two elements in 'arrays' object
				stdevs[i] = Math.sqrt(stdevs[i]);
			}			
		}
		return stdevs;		
	}/*  end of calculateStandardDeviations method  */
	
	
	/*  start of outputArrayToFile method  */
	private void outputArrayToFile(double[] array, String outputFilename){
		try{
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + outputFilename));
			
			if(array.length == 0){
				bw.write("Note: the array to output had size zero (0). There is no data to output for this array.");
				bw.flush();
			}
			else{
				for(int i=0; i<array.length; i=i+resolution){
					bw.write((i+minBase) + "\t" + array[i]);
					bw.newLine();
					bw.flush();
				}
			}
			bw.flush();
			bw.close();
			System.out.println("\tOutput file: " + outputPath + outputFilename);			
		}catch (IOException e) {
			e.printStackTrace();
		}//end of try-catch block
	}/*  end of outputArrayToFile method  */
	
	
	/*  start of calculateBaseDensityStandAlone method  */
	private void calculateBaseDensityStandAlone(){
		baseDensities = new ArrayList<double[]>(numberOfStructures);
		
		
		//  calculate the base density for each of the structures individually
		for(int structureIndex=0; structureIndex<numberOfStructures; structureIndex++){
			System.out.println("\tWorking on structure " + (structureIndex+1) + " of " + numberOfStructures);
			double[] currentStructureBaseDensities = new double[numberOfBases];
			
			//  go through every pair of bases in the structure
			for(int firstBaseIndex=minBase; firstBaseIndex<maxBase+1-resolution; firstBaseIndex=firstBaseIndex+resolution){
				for(int secondBaseIndex=firstBaseIndex+resolution; secondBaseIndex<maxBase+1; secondBaseIndex=secondBaseIndex+resolution){
					//  update the base density for each density base that is contained within the radius of the bounding sphere
					if(getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) < radiusSquared){					
						currentStructureBaseDensities[firstBaseIndex-minBase]++;
						currentStructureBaseDensities[secondBaseIndex-minBase]++;
					}				
				}
			}
			baseDensities.add(currentStructureBaseDensities);
		}
		baseDensities.trimToSize();
		

		//  calculate the average base density and standard deviation
		averageBaseDensity = calculateAverages(baseDensities);
		stdevBaseDensity = calculateStandardDeviations(baseDensities, averageBaseDensity);
		
		
		//  output the average base density and standard deviation to file
		outputArrayToFile(averageBaseDensity, "averageBaseDensity.txt");
		outputArrayToFile(stdevBaseDensity, "stdevBaseDensity.txt");
	}/*  end of calculateBaseDensityStandAlone method  */
	
	
	/*  start of calculateBaseCondensationStandAlone method  */
	private void calculateBaseCondensationStandAlone(){
		baseCondensations = new ArrayList<double[]>(numberOfStructures);
		
		
		//  calculate the base condensation for each of the structures individually
		for(int structureIndex=0; structureIndex<numberOfStructures; structureIndex++){
			System.out.println("\tWorking on structure " + (structureIndex+1) + " of " + numberOfStructures);
			
			double[] currentStructureBaseCondensations = new double[numberOfBases];
			
			//  go through every pair of bases in the structure
			for(int firstBaseIndex=minBase; firstBaseIndex<maxBase+1-resolution; firstBaseIndex=firstBaseIndex+resolution){
				for(int secondBaseIndex=firstBaseIndex+resolution; secondBaseIndex<maxBase+1; secondBaseIndex=secondBaseIndex+resolution){
					//  update the base condensation for each condensation base that is contained within the radius of the bounding sphere
					if(getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) > radiusSquared){
						break;
					}
					else{
						currentStructureBaseCondensations[firstBaseIndex-minBase]++;
						currentStructureBaseCondensations[secondBaseIndex-minBase]++;
					}
				}
			}
			baseCondensations.add(currentStructureBaseCondensations);
		}
		baseCondensations.trimToSize();
		
		
		//  calculate the average base condensation and standard deviation
		averageBaseCondensation = calculateAverages(baseCondensations);
		stdevBaseCondensation = calculateStandardDeviations(baseCondensations, averageBaseCondensation);
		
		
		//  output the average base condensation and standard deviation to file
		outputArrayToFile(averageBaseCondensation, "averageBaseCondensation.txt");
		outputArrayToFile(stdevBaseCondensation, "stdevBaseCondensation.txt");
	}/*  end of calculateBaseCondensationStandAlone method  */
	
	
	/*  start of calculateBaseLoopingStandAlone method  */
	private void calculateBaseLoopingStandAlone(){
		baseLoopings = new ArrayList<double[]>(numberOfStructures);
		
		
		//  calculate the base looping for each of the structures individually
		for(int structureIndex=0; structureIndex<numberOfStructures; structureIndex++){
			System.out.println("\tWorking on structure " + (structureIndex+1) + " of " + numberOfStructures);
			double[] currentStructureBaseLooping = new double[numberOfBases];
			
			//  go through every pair of bases in the structure
			for(int firstBaseIndex=minBase; firstBaseIndex<maxBase+1-resolution; firstBaseIndex=firstBaseIndex+resolution){
				boolean condensationPhaseOver = false;
				for(int secondBaseIndex=firstBaseIndex+resolution; secondBaseIndex<maxBase+1; secondBaseIndex=secondBaseIndex+resolution){
					if(!condensationPhaseOver && (getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) > radiusSquared)){
						condensationPhaseOver = true;
					}
					//  update the base looping for each looping base that is contained within the radius of the bounding sphere
					if(condensationPhaseOver && (getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) < radiusSquared)){					
						currentStructureBaseLooping[firstBaseIndex-minBase]++;
						currentStructureBaseLooping[secondBaseIndex-minBase]++;
					}				
				}
			}
			baseLoopings.add(currentStructureBaseLooping);
		}
		baseLoopings.trimToSize();
		

		//  calculate the average base looping and standard deviation
		averageBaseLooping = calculateAverages(baseLoopings);
		stdevBaseLooping = calculateStandardDeviations(baseLoopings, averageBaseLooping);
		
		
		//  output the average base looping and standard deviation to file
		outputArrayToFile(averageBaseLooping, "averageBaseLooping.txt");
		outputArrayToFile(stdevBaseLooping, "stdevBaseLooping.txt");
	}/*  end of calculateBaseLoopingStandAlone method  */
	

	/*  start of calculateBaseDensityAsSumOfCondensationAndLooping method  */
	private void calculateBaseDensityAsSumOfCondensationAndLooping(){
		baseDensities = new ArrayList<double[]>(numberOfStructures);
		
		
		//  calculate the base density for each of the structures individually
		for(int structureIndex=0; structureIndex<numberOfStructures; structureIndex++){
			System.out.println("\tWorking on structure " + (structureIndex+1) + " of " + numberOfStructures);
			double[] currentStructureBaseDensities = new double[numberOfBases];
			
			//  go through every base in the structure
			for(int i=0; i<numberOfBases; i=i+resolution){
				//  update the base density for each base					
				currentStructureBaseDensities[i] = baseCondensations.get(structureIndex)[i] + baseLoopings.get(structureIndex)[i];
			}
			baseDensities.add(currentStructureBaseDensities);
		}
		baseDensities.trimToSize();
		

		//  calculate the average base density and standard deviation
		averageBaseDensity = calculateAverages(baseDensities);
		stdevBaseDensity = calculateStandardDeviations(baseDensities, averageBaseDensity);
		
		
		//  output the average base density and standard deviation to file
		outputArrayToFile(averageBaseDensity, "averageBaseDensity.txt");
		outputArrayToFile(stdevBaseDensity, "stdevBaseDensity.txt");
	}/*  end of calculateBaseDensityAsSumOfCondensationAndLooping method  */
	
	
	/*  start of threeMusketeers method  */
	private void threeMusketeers(){
		baseCondensations = new ArrayList<double[]>(numberOfStructures);
		baseLoopings = new ArrayList<double[]>(numberOfStructures);
		baseDensities = new ArrayList<double[]>(numberOfStructures);
		
		//  calculate the base condensation, looping, and density for each of the structures individually
		for(int structureIndex=0; structureIndex<numberOfStructures; structureIndex++){
			System.out.println("\tWorking on structure " + (structureIndex+1) + " of " + numberOfStructures);
			
			double[] currentStructureBaseCondensations = new double[numberOfBases];
			double[] currentStructureBaseLooping = new double[numberOfBases];
			double[] currentStructureBaseDensities = new double[numberOfBases];
			
			//  go through every pair of bases in the structure
			for(int firstBaseIndex=minBase; firstBaseIndex<maxBase+1-resolution; firstBaseIndex=firstBaseIndex+resolution){
				boolean condensationPhaseOver = false;
				for(int secondBaseIndex=firstBaseIndex+resolution; secondBaseIndex<maxBase+1; secondBaseIndex=secondBaseIndex+resolution){
					//  update the base condensation and base looping for each base
					
					if(!condensationPhaseOver){
						if(getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) > radiusSquared){
							condensationPhaseOver = true;
						}
						else{
							//  update the base condensation for each base that is contained within the radius of the bounding sphere
							currentStructureBaseCondensations[firstBaseIndex-minBase]++;
							currentStructureBaseCondensations[secondBaseIndex-minBase]++;
						}
					}
					//  update the base looping for each base that is contained within the radius of the bounding sphere
					else if(getDistanceSquaredBetweenBases(structureIndex, firstBaseIndex, secondBaseIndex) < radiusSquared){					
						currentStructureBaseLooping[firstBaseIndex-minBase]++;
						currentStructureBaseLooping[secondBaseIndex-minBase]++;
					}
				}
			}
			//  update the base density using the sum of the base condensation and base looping method
			for(int i=0; i<numberOfBases; i=i+resolution){
				currentStructureBaseDensities[i] = currentStructureBaseCondensations[i] + currentStructureBaseLooping[i];
			}
			baseCondensations.add(currentStructureBaseCondensations);
			baseLoopings.add(currentStructureBaseLooping);
			baseDensities.add(currentStructureBaseDensities);
		}
		baseCondensations.trimToSize();
		baseLoopings.trimToSize();
		baseDensities.trimToSize();
		
		
		//  calculate the averages and standard deviations
		averageBaseCondensation = calculateAverages(baseCondensations);
		stdevBaseCondensation = calculateStandardDeviations(baseCondensations, averageBaseCondensation);
		averageBaseLooping = calculateAverages(baseLoopings);
		stdevBaseLooping = calculateStandardDeviations(baseLoopings, averageBaseLooping);
		averageBaseDensity = calculateAverages(baseDensities);
		stdevBaseDensity = calculateStandardDeviations(baseDensities, averageBaseDensity);
		
		
		//  output the averages and standard deviations to file
		outputArrayToFile(averageBaseCondensation, "averageBaseCondensation.txt");
		outputArrayToFile(stdevBaseCondensation, "stdevBaseCondensation.txt");
		outputArrayToFile(averageBaseLooping, "averageBaseLooping.txt");
		outputArrayToFile(stdevBaseLooping, "stdevBaseLooping.txt");
		outputArrayToFile(averageBaseDensity, "averageBaseDensity.txt");
		outputArrayToFile(stdevBaseDensity, "stdevBaseDensity.txt");
	}/*  end of threeMusketeers method  */
		
	
	/*  start of run method  */
	public void run(String[] args){
		if(args.length == 5){
			experimentName = args[0];
			fragmentsName = args[1];
			numberOfStructures = Integer.parseInt(args[2]);
			resolution = Integer.parseInt(args[3]);
			radius = Double.parseDouble(args[4]);
		}
		else{
			System.out.println("Error:\tplease use the following pattern (without the quotes)\n" +
					"\t\"java MicrocosmMain experimentName fragmentsFilename.txt numberOfStructures resolution radius\"\n" +
					"Notes:\tthe default numberOfStructures is " + numberOfStructures + "\n" +
					"\tthe default resolution is " + resolution + "\n" +
					"\tthe default radius is " + radius + "\n" +
					"\tNote: if there is no experiment name then specify \"none\" as the experiment name\n" +
					"System.exit(0)");
			System.exit(0);
		}
		
		
		//  create the output folder
		if(!experimentName.equals("none")){
			outputFolder += "_" + experimentName;
		}
		outputPath = outputFolder + "/";
		new File(outputFolder).mkdir();
		
		
		//  display the parameters of the run
		System.out.println("Parameters:");
		if(numberOfStructures == 1){
			System.out.println("\tUsing a pool of " + numberOfStructures + " structure");
		}
		else{
			System.out.println("\tUsing a pool of " + numberOfStructures + " structures");
		}
		System.out.println("\tresolution = " + resolution);
		System.out.println("\tradius = " + radius);
		
		
		//  output the parameters to file
		try{
			String parameterFilename = "parameters.txt";
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + parameterFilename));
			
			bw.write("experimentName=\t" + experimentName);
			bw.newLine();
			bw.flush();
			bw.write("fragmentsName=\t" + fragmentsName);
			bw.newLine();
			bw.flush();
			bw.write("numberOfStructures=\t" + numberOfStructures);
			bw.newLine();
			bw.flush();
			bw.write("resolution=\t" + resolution);
			bw.newLine();
			bw.flush();
			bw.write("radius=\t" + radius);
			bw.newLine();
			bw.flush();
			bw.write("command line parameters:");
			for(String arg: args){
				bw.write("\t" + arg);
			}
			bw.newLine();
			bw.flush();
			
			bw.close();
			System.out.println("\tOutput file: " + outputPath + parameterFilename);
		}catch (IOException e) {
			e.printStackTrace();
		}//end of try-catch block
		
		
		System.out.println("\nRunning...\n");
		
		
		//  set the radius squared value
		radiusSquared = radius*radius;
		
		
		//  calculate how long running Microcosm takes (time point before run)
		long timeBefore = System.currentTimeMillis();

		
		//  get the number of points in each structure
		//  Note: need to get the number of points before reading in the fragments file
		if(experimentName.equals("none")){
			numberOfPoints = pdb.getNumberOfPoints(inputPath + "structure_0001" + ".pdb");
		}
		else{
			numberOfPoints = pdb.getNumberOfPoints(inputPath + "structure_0001" + "_" + experimentName + ".pdb");
		}
//		System.out.println("\tNumber of points = " + numberOfPoints);
		
		
		//  load all of the structures from file
		getStructuresFromFile();
//		System.out.println("\tstructures.size() = " + structures.size());
		
		
		//  load all of the fragments from file
		getFragmentsFromFile();
		
		
		//  calculate base density using the stand alone method
//		System.out.println("Calculating the base density (stand alone method)...");
//		calculateBaseDensityStandAlone();
//		System.out.println("Done calculating the base density!\n");
		
		
		//  calculate base condensation using the stand alone method
//		System.out.println("Calculating the base condensation (stand alone method)...");
//		calculateBaseCondensationStandAlone();
//		System.out.println("Done calculating the base condensation!\n");
		
		
		//  calculate base looping using the stand alone method
//		System.out.println("Calculating the base looping (stand alone method)...");
//		calculateBaseLoopingStandAlone();
//		System.out.println("Done calculating the base looping!\n");
		
		
		//  calculate base density using the sum of the condensation and looping method
		//  Note: this method can only be called after both the base condensation and base looping methods have been called
//		System.out.println("Calculating the base density (sum of condensation and looping method)...");
//		calculateBaseDensityAsSumOfCondensationAndLooping();
//		System.out.println("Done calculating the base density!\n");

		
		//  calculate base condensation, base looping, and base density using the three musketeers method
		System.out.println("Calculating the base condensation, base looping, and base density (three musketeers method)...");
		threeMusketeers();
		System.out.println("Done calculating the base condensation, base looping, and base density!\n");
		
		
		//  output the time taken by Microcosm
		long timeTaken = System.currentTimeMillis() - timeBefore;
		timeTaken = timeTaken / (long)(1000.0);			//	convert from milliseconds to seconds		
		if(timeTaken > 3600){		//	(1hr = 3600 sec)
			int numberOfHours = (int)timeTaken / (int)3600;
			long leftOver = timeTaken % 3600;
			int numberOfMinutes = (int)leftOver / (int)60;
			leftOver = leftOver % 60;
			int numberOfSeconds = (int)leftOver;
			System.out.println("Amount of time taken = " + numberOfHours + "hr" + numberOfMinutes + "min" + numberOfSeconds + "sec");
		}
		else if(timeTaken > 60){		//	(1min = 60 sec)
			int numberOfMinutes = (int)timeTaken / (int)60;
			long leftOver = timeTaken % 60;
			int numberOfSeconds = (int)leftOver;
			System.out.println("Amount of time taken = " + numberOfMinutes + "min" + numberOfSeconds + "sec");
		}
		else{
			int numberOfSeconds = (int)timeTaken;
			System.out.println("Amount of time taken = " + numberOfSeconds + "sec");
		}
		
		
		System.out.print("\nDone!");
	}/*  end of run method  */
}/*  end of Microcosm class  */
